function corr_XY = Calc_quarterly_corr(Phi_Nz,stationary_prob_Nz,X_Nz,Y_Nz)
% computes:
% corr((X(t+1)+X(t+2)+X(t+3))/3,(Y(t+1)+Y(t+2)+Y(t+3))/3)
% E[(X(t+1)+X(t+2)+X(t+3))/3]
% Var[(X(t+1)+X(t+2)+X(t+3))/3]

    % E[X], E[Y]s
    [X_mean, X_std] = Calc_quarterly_ave_NHz(Phi_Nz,stationary_prob_Nz,X_Nz);
    [Y_mean, Y_std] = Calc_quarterly_ave_NHz(Phi_Nz,stationary_prob_Nz,Y_Nz);

    % compute cross terms
    XY_cross_terms = sum(stationary_prob_Nz(:).*(X_Nz(:).*( 3*Y_Nz(:) + Phi_Nz*(2*Y_Nz(:) + Phi_Nz*Y_Nz(:)) ) ...
        + Y_Nz(:).*( Phi_Nz*(2*X_Nz(:) + Phi_Nz*X_Nz(:)) ) ));

    corr_XY = (XY_cross_terms/9 - X_mean*Y_mean)/(X_std*Y_std);

end